Bayesian statistical modeling by Stan&R introduction (the book by baba)
01
data was downloaded from this website
library(tidyverse)
library(rstan)
library(bayesplot)
path <- "/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-4-1-beer-sales-1.csv"
df <- read_csv(path)
df# A tibble: 100 x 1
sales
<dbl>
1 87.5
2 104.
3 83.3
4 132.
5 107.
6 83.6
7 110.
8 115.
9 112.
10 93.9
# … with 90 more rowsmcmc_result <- stan(file = "var-mean.stan", data = list(sales = df$sales, N = nrow(df)),
seed = 1, chains = 4, iter = 2000, warmup = 1000, thin = 1)
SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.3 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.053426 seconds (Warm-up)
Chain 1: 0.03644 seconds (Sampling)
Chain 1: 0.089866 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.03837 seconds (Warm-up)
Chain 2: 0.035696 seconds (Sampling)
Chain 2: 0.074066 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 8e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.048032 seconds (Warm-up)
Chain 3: 0.031268 seconds (Sampling)
Chain 3: 0.0793 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'var-mean' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.048851 seconds (Warm-up)
Chain 4: 0.064571 seconds (Sampling)
Chain 4: 0.113422 seconds (Total)
Chain 4: mcmc_sample <- rstan::extract(mcmc_result, permuted = FALSE)
mcmc_sample, , parameters = mu
chains
iterations chain:1 chain:2 chain:3 chain:4
[1,] 100.64870 104.79868 100.68753 101.61185
[2,] 102.79747 105.97725 102.40376 101.37497
[3,] 100.64631 101.67638 101.53096 102.84214
[4,] 100.56696 101.09195 102.40958 103.42633
[5,] 102.79871 103.24347 104.02485 103.08640
[6,] 102.69393 103.24347 99.64743 100.91259
[7,] 101.57675 102.21287 97.71370 98.74514
[8,] 102.89931 100.23452 103.89675 104.65942
[9,] 99.62106 99.40103 101.03884 102.73738
[10,] 103.55490 99.40626 102.43715 102.94000
[11,] 102.90797 103.55604 102.74112 102.62739
[12,] 103.86359 99.09045 104.93339 102.48959
[13,] 101.24143 99.09045 100.14706 104.38543
[14,] 99.85707 105.49999 99.87357 99.39213
[15,] 101.37068 105.18203 100.20999 105.90998
[16,] 103.57256 104.52472 101.36059 103.34531
[17,] 101.92233 102.00397 105.07658 100.00039
[18,] 103.01964 101.83571 101.04252 103.92699
[ reached getOption("max.print") -- omitted 982 row(s) and 2 matrix slice(s) ]mcmc_sample[1, "chain:3", "mu"][1] 100.6875mcmc_sample[, "chain:1", "mu"] [1] 100.64870 102.79747 100.64631 100.56696 102.79871 102.69393 101.57675
[8] 102.89931 99.62106 103.55490 102.90797 103.86359 101.24143 99.85707
[15] 101.37068 103.57256 101.92233 103.01964 104.38627 101.61651 103.18445
[22] 100.15041 104.19408 101.60059 100.02435 103.37332 104.80943 101.02734
[29] 99.75564 100.35635 100.97986 101.48021 101.62266 101.94145 101.89881
[36] 102.81174 101.82513 102.36507 101.81274 102.49252 101.28915 99.51839
[43] 99.35342 98.70641 102.85533 101.80770 101.55767 102.99445 101.20891
[50] 102.10344 104.19419 104.39294 101.11485 99.07752 100.25881 100.33359
[57] 100.01780 100.70371 99.51196 100.53404 103.08729 102.13608 103.13025
[64] 106.22945 98.28900 104.21365 103.10044 101.40143 103.48623 103.35901
[71] 104.54869 104.36264 103.02141 102.05826 101.08636
[ reached getOption("max.print") -- omitted 925 entries ]mcmc_sample[, , "mu"] chains
iterations chain:1 chain:2 chain:3 chain:4
[1,] 100.64870 104.79868 100.68753 101.61185
[2,] 102.79747 105.97725 102.40376 101.37497
[3,] 100.64631 101.67638 101.53096 102.84214
[4,] 100.56696 101.09195 102.40958 103.42633
[5,] 102.79871 103.24347 104.02485 103.08640
[6,] 102.69393 103.24347 99.64743 100.91259
[7,] 101.57675 102.21287 97.71370 98.74514
[8,] 102.89931 100.23452 103.89675 104.65942
[9,] 99.62106 99.40103 101.03884 102.73738
[10,] 103.55490 99.40626 102.43715 102.94000
[11,] 102.90797 103.55604 102.74112 102.62739
[12,] 103.86359 99.09045 104.93339 102.48959
[13,] 101.24143 99.09045 100.14706 104.38543
[14,] 99.85707 105.49999 99.87357 99.39213
[15,] 101.37068 105.18203 100.20999 105.90998
[16,] 103.57256 104.52472 101.36059 103.34531
[17,] 101.92233 102.00397 105.07658 100.00039
[18,] 103.01964 101.83571 101.04252 103.92699
[ reached getOption("max.print") -- 982 行を無視しました ] mu_mcmc_vec <- as.vector(mcmc_sample[, , "mu"])
mu_mcmc_vec [1] 100.64870 102.79747 100.64631 100.56696 102.79871 102.69393 101.57675
[8] 102.89931 99.62106 103.55490 102.90797 103.86359 101.24143 99.85707
[15] 101.37068 103.57256 101.92233 103.01964 104.38627 101.61651 103.18445
[22] 100.15041 104.19408 101.60059 100.02435 103.37332 104.80943 101.02734
[29] 99.75564 100.35635 100.97986 101.48021 101.62266 101.94145 101.89881
[36] 102.81174 101.82513 102.36507 101.81274 102.49252 101.28915 99.51839
[43] 99.35342 98.70641 102.85533 101.80770 101.55767 102.99445 101.20891
[50] 102.10344 104.19419 104.39294 101.11485 99.07752 100.25881 100.33359
[57] 100.01780 100.70371 99.51196 100.53404 103.08729 102.13608 103.13025
[64] 106.22945 98.28900 104.21365 103.10044 101.40143 103.48623 103.35901
[71] 104.54869 104.36264 103.02141 102.05826 101.08636
[ reached getOption("max.print") -- omitted 3925 entries ]mean(mu_mcmc_vec)[1] 102.166quantile(mu_mcmc_vec, probs = c(0.025, 0.975)) 2.5% 97.5%
98.54126 105.75601 print(mcmc_result, probs = c(0.025, 0.5, 0.975))Inference for Stan model: var-mean.
4 chains, each with iter=2000; warmup=1000; thin=1;
post-warmup draws per chain=1000, total post-warmup draws=4000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
mu 102.17 0.03 1.82 98.54 102.19 105.76 2889 1
sigma 18.18 0.02 1.30 15.90 18.09 20.98 3296 1
lp__ -336.44 0.02 0.97 -339.06 -336.15 -335.47 1877 1
Samples were drawn using NUTS(diag_e) at Thu Jan 2 19:00:40 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).library(ggfortify)
autoplot(ts(mcmc_sample[, , "mu"]), facets = F, ylab = "mu", main = "trace plot") +
theme_bw(base_size = 14)
ts(mcmc_sample[, , "mu"])Time Series:
Start = 1
End = 1000
Frequency = 1
chain:1 chain:2 chain:3 chain:4
1 100.64870 104.79868 100.68753 101.61185
2 102.79747 105.97725 102.40376 101.37497
3 100.64631 101.67638 101.53096 102.84214
4 100.56696 101.09195 102.40958 103.42633
5 102.79871 103.24347 104.02485 103.08640
6 102.69393 103.24347 99.64743 100.91259
7 101.57675 102.21287 97.71370 98.74514
8 102.89931 100.23452 103.89675 104.65942
9 99.62106 99.40103 101.03884 102.73738
10 103.55490 99.40626 102.43715 102.94000
11 102.90797 103.55604 102.74112 102.62739
12 103.86359 99.09045 104.93339 102.48959
13 101.24143 99.09045 100.14706 104.38543
14 99.85707 105.49999 99.87357 99.39213
15 101.37068 105.18203 100.20999 105.90998
16 103.57256 104.52472 101.36059 103.34531
17 101.92233 102.00397 105.07658 100.00039
18 103.01964 101.83571 101.04252 103.92699
[ reached getOption("max.print") -- 982 行を無視しました ] ggplot(data = mcmc_sample[, , "mu"] %>% as.vector %>% data.frame(sample = .)) +
geom_density(aes(x = sample), fill = "skyblue3", color = "skyblue4", alpha = 0.6) +
theme_bw(base_size = 14)
mcmc_hist(mcmc_sample, pars = c("mu", "sigma"))
mcmc_dens(mcmc_sample, pars = c("mu", "sigma"))
mcmc_trace(mcmc_sample, pars = c("mu", "sigma"))
mcmc_combo(mcmc_sample, pars = c("mu", "sigma"))
mcmc_intervals(mcmc_sample, pars = c("mu", "sigma"), prob = 0.8, prob_outer = 0.95)
mcmc_areas(mcmc_sample, pars = c("mu", "sigma"), prob = 0.8, prob_outer = 0.95)
mcmc_acf_bar(mcmc_sample, pars = c("mu", "sigma"))
animal_num <- read_csv("/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-5-1-animal-num.csv")
animal_num %>% table.
0 1 2 3 4
64 85 39 8 4 mcmc_normal <- stan(file = "animal_normal_dist.stan", data = list(animal_num = animal_num$animal_num,
N = nrow(animal_num)), seed = 1, chains = 4, iter = 2000, warmup = 1000,
thin = 1)
SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.4e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.04218 seconds (Warm-up)
Chain 1: 0.033597 seconds (Sampling)
Chain 1: 0.075777 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 1.2e-05 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.05225 seconds (Warm-up)
Chain 2: 0.033195 seconds (Sampling)
Chain 2: 0.085445 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 8e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.038179 seconds (Warm-up)
Chain 3: 0.04208 seconds (Sampling)
Chain 3: 0.080259 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'animal_normal_dist' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.035994 seconds (Warm-up)
Chain 4: 0.035385 seconds (Sampling)
Chain 4: 0.071379 seconds (Total)
Chain 4: mcmc_poisson <- stan(file = "animal_poisson_dist.stan", data = list(animal_num = animal_num$animal_num,
N = nrow(animal_num)), seed = 1, chains = 4, iter = 2000, warmup = 1000,
thin = 1)
SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1.3e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.035789 seconds (Warm-up)
Chain 1: 0.034036 seconds (Sampling)
Chain 1: 0.069825 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 9e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.033882 seconds (Warm-up)
Chain 2: 0.028672 seconds (Sampling)
Chain 2: 0.062554 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 6e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.06 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.031612 seconds (Warm-up)
Chain 3: 0.027575 seconds (Sampling)
Chain 3: 0.059187 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'animal_poisson_dist' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 8e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.08 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.028679 seconds (Warm-up)
Chain 4: 0.026872 seconds (Sampling)
Chain 4: 0.055551 seconds (Total)
Chain 4: - prediction
y_rep_normal <- extract(mcmc_normal)$pred
y_rep_poisson <- extract(mcmc_poisson)$preddim(y_rep_normal)[1] 4000 200y_rep_normal[1, ] [1] 1.146980013 1.854437914 1.021145357 0.578413718 0.258020066
[6] -0.351310680 1.253304848 0.110552641 0.755742723 1.228747160
[11] 1.344334788 0.541145278 0.780954016 1.949987041 0.043828549
[16] 0.721685533 2.040762643 1.010444509 1.047763208 2.621768586
[21] 1.742483477 1.983092709 2.512677485 1.647432538 -0.827800346
[26] 2.267476210 0.883656061 0.629227036 0.947072106 3.769873426
[31] -0.003849530 0.688832818 1.612792468 1.174575527 1.192150488
[36] 2.077392144 1.018682377 0.629081869 2.157806212 1.411999471
[41] 1.154246517 1.852153406 0.779993220 1.153493484 -0.686172286
[46] 1.687172205 1.037327786 0.029628722 1.555700961 1.719985511
[51] 1.215092548 -0.002299914 0.521871637 1.369611267 0.558290493
[56] -0.855950558 1.335117234 2.746824217 0.752569171 1.938582522
[61] 0.588832991 0.801822983 3.518490204 -0.113727357 1.036300359
[66] 0.718570796 0.640913945 -0.140313360 0.621811302 2.642044739
[71] 0.823595780 0.658087706 0.454968067 0.845824208 1.100157136
[ reached getOption("max.print") -- omitted 125 entries ]y_rep_poisson[1, ] [1] 1 1 2 1 0 0 1 0 1 2 3 0 1 3 1 1 1 0 0 1 1 2 3 6 1 0 1 5 0 0 0 3 0 3 0
[36] 0 4 0 1 1 0 0 2 2 1 1 0 1 1 1 1 0 2 0 0 0 1 0 1 0 1 1 0 0 0 0 3 1 1 0
[71] 2 1 0 1 1
[ reached getOption("max.print") -- omitted 125 entries ]ppc_hist(y = animal_num$animal_num, yrep = y_rep_normal[1:5, ])
path <- "/Users/tomoya/Documents/dataset/bayesian_intro_baba/2-6-1-beer-sales-ab.csv"
sales_ab <- read_csv(path)ggplot(data = sales_ab, aes(x = sales, y = ..density.., color = beer_name, fill = beer_name)) +
geom_histogram(alpha = 0.5, position = "identity") + geom_density(alpha = 0.5,
size = 0) + theme_bw(base_size = 14)
inspect <- sales_ab %>% mosaic::inspect()
inspect$categorical# A tibble: 1 x 6
name class levels n missing distribution
<chr> <chr> <int> <int> <int> <chr>
1 beer_na… charact… 2 200 0 "A (50%), B (50%) …inspect$quantitative# A tibble: 1 x 11
name class min Q1 median Q3 max mean sd n missing
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <int>
1 sales numeric 55.7 103. 128. 165. 239. 136. 41.1 200 0data_list_ab <- list(sales_a = sales_ab %>% filter(beer_name == "A") %>% select(sales) %>%
1[[]], sales_b = sales_ab %>% filter(beer_name == "B") %>% select(sales) %>%
1[[]], N = 100)mcmc_result_sales_ab <- stan(file = "sales_ab_normal.stan", data = data_list_ab,
seed = 1, chains = 4, iter = 2000, warmup = 1000, thin = 1)
SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 1e-05 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.1 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 0.050558 seconds (Warm-up)
Chain 1: 0.019665 seconds (Sampling)
Chain 1: 0.070223 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 2).
Chain 2:
Chain 2: Gradient evaluation took 5e-06 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2:
Chain 2:
Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 2:
Chain 2: Elapsed Time: 0.040216 seconds (Warm-up)
Chain 2: 0.021298 seconds (Sampling)
Chain 2: 0.061514 seconds (Total)
Chain 2:
SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 3).
Chain 3:
Chain 3: Gradient evaluation took 5e-06 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.05 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3:
Chain 3:
Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 3:
Chain 3: Elapsed Time: 0.037421 seconds (Warm-up)
Chain 3: 0.04389 seconds (Sampling)
Chain 3: 0.081311 seconds (Total)
Chain 3:
SAMPLING FOR MODEL 'sales_ab_normal' NOW (CHAIN 4).
Chain 4:
Chain 4: Gradient evaluation took 4e-06 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.04 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4:
Chain 4:
Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
Chain 4:
Chain 4: Elapsed Time: 0.047191 seconds (Warm-up)
Chain 4: 0.04551 seconds (Sampling)
Chain 4: 0.092701 seconds (Total)
Chain 4: mcmc_sample <- extract(mcmc_result_sales_ab, permuted = F)
mcmc_dens(mcmc_sample, pars = "diff") + theme_bw(base_size = 14)